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Abstract 

We describe a rapidly converging algorithm for solving the Kohn-Sham equations 
and equations of similar structure that appear frequently in calculations of the 
structure of inhomogeneous electronic many-body systems. The algorithm has its 
roots the Hohenberg-Kohn theorem and solves directly for the electron density; 
single-particle wave functions are only used as auxiliary quantities. The method 
has been implemented for symmetric "slabs" of jellium as well as for spherical 
jellium clusters. Starting from very rough guesses for the initial electron density, 
convergence is reached within a few iterations. The iterations are driven by the 
static electric susceptibility. 



1 INTRODUCTION 



Density-functional theory [1] is a popular method for studying properties of 
inhomogeneous many-electron systems. Central to the execution of the theory 
is the solution of the Kohn-Sham equation, which is an effective Schrodinger 
equation for a set of single-particle wave functions 0i(r) 

- |^V 2 0i(r) + VW](r)0,(r) = e^(r) , (1.1) 



which are characterized by a set of quantum numbers i. Vh (r) in Eq. (1.1) is an 
effective one-body potential that describes both the influence of an external 
field and many-body interactions and depends implicitly on the density of 
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the system. In density functional theory, this effective one-body ( "Hartree" )- 
potential is 



V h (t)=V c (t)+V xc {t) 



(1.2) 



where Vc(t) is the Coulomb potential 



Vc(r) = fd 3 r'-^-- [p(r') - p + (r')} , 



(1.3) 



and V xc (r) is the "exchange-correlation" correction. p+(r') is the charge den- 
sity of the positive background. 

From the single-particle wave functions 0i(r) one constructs the physical one- 
body density 



the n(i) are the occupation numbers of the single-particle states. 

Structurally similar equations also appear in microscopic theories of inhomo- 
geneous, correlated electrons [2,3]. In fact, equations equivalent to the Kohn- 
Sham equation appear in the theory of basically all inhomogeneous quantum 
many-particle systems such as nuclei, quantum liquid clusters, and quantum 
liquid films. We disregard, for the time being, non-localities like a Fock-term 
and will comment on how to include such complications further below. 

In this paper we will develop a rapidly converging iterative scheme that solves, 
very much in the spirit of the Hohenberg-Kohn theorem, directly for the den- 
sity p(r). 

It is plausible to try to solve Eqs. (1.1)-(1.4) iteratively: One calculates, from 
an initial guess for the density p(r), the single-particle orbitals 0i(r) from Eq. 
(1.1), from these a new single-particle density (1.4), and iterates the proce- 
dure until convergence is reached. This procedure converges well for atoms 
and other systems that are dominated by an external field, but the conver- 
gence can be very slow for metallic clusters, self-bound systems, or in cases 
where local charge-neutrality is important. Judicious "mixing" of consecutive 
iterations is needed [4,5] with admixtures of the "new" solutions of less than 
1 percent for clusters of a few hundred electrons. A specialized "gradient iter- 
ation method" [6,7] improves the stability of the process considerably, but 
still takes very many iterations to converge. This situation is unsatisfactory 
because one should think that, because these systems are very stable in na- 



p(r) = 5>(i)|0i(r)| 2 ; 



(1.4) 
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ture, one should be able to mimic nature and to design a stable and rapidly 
converging numerical algorithm to obtain ground state configurations. 

The cause for the sometimes delicate convergence properties is rather clear in 
electronic systems: The configuration with local charge neutrality will mini- 
mize the potential energy which can, on the other hand, be enormous if local 
charge neutrality is not maintained. The quantum kinetic energy causes only 
a small violation of local charge neutrality. In other words, the quantity that 
drives the configuration towards the correct ground state is the potential en- 
ergy, whereas the simple iterative procedure (1.1), (1.4) focuses on the kinetic 
energy. 



2 Algorithm 

Our procedure to solve the coupled equation Eq. (1.1)— (1.4) has its roots in 
the Hohenberg-Kohn theorem [8] that states, among others, that the one- 
body density p(r) is the only truly independent variable. We therefore design 
an algorithm that solves, by a Newton- Raphson procedure, directly for p(r). 

Consider the functional 



where the 0i[p](r) are the solutions of Eq. (1.1) for the density p(r). The exact 
solution of Eq. (1.1) is the density that satisfies 



Non-linear equations of the type (2.2) can be solved iteratively by the Newton- 
Raphson algorithm. Assume that we have an k-th guess for the density, say 
p( fc )(r). A next iteration p( fc+1 )(r) = p^(r)+Sp^ k \r) is defined by the condition 



F[p](r) = £Mi)l0i[p](r)| 2 -p(r) 



(2.1) 



F[p](r) = 0. 



(2.2) 



= F[p( fc+1 )](r) = F[pW + 5p^](r) 



= F[pW](r) + J d 3 r' 




5p^(r')+0((5p^) 2 ) . 



(2.3) 



Defining 



SF[p](r) 
Sp(r f ) 



e(r,r';0), 



(2.4) 
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the density correction 5p^(r) is determined by 



F[pW](r) = J d 3 r'e(r, r'; 0)5p (fc) (r) • (2-5) 

which is, upon discretization, a system of linear equations which can be solved 
by standard methods. 

The remaining task is the calculation and physical interpretation of the kernel 
e(r, r';0). For any set of single particle orbitals {0i(r)} representing a k-ih 
iteration {(j)[ k \r)}, we have 

e(r,r';0)=5(r-r')-£n(i) 



0i« 



8p{r>) 



+ c. c. 



(2.6) 



The variations of the single-particle wave functions are obtained by first order 
perturbation theory: 

<*fr(r) _ f « r n y WjjO ^(Q r2? , 
(Jp(r') " 7 ^ £J - e, 5p(r') " 1 " > 



Inserting (2.7) in the second term of Eq. (2.6), we find 



Vr" 



>: r 



5p(r') 



+ c. c. 



2^n(i)n(j) =! + c. c. 



e, - £i 



W g (r") 
5p(r') 



J dV'xo(r,r"; 0)F p _ h (r", r') ■ 



(2.8) 



Two new quantities have been introduced in the last line of Eq. (2.8): In 
the density variation of the Hartree-potential we recover a local and static 
approximation for the effective "particle-hole" interaction, 



Vh(r,r') = 



SVh(t) 
5p(r>) 



(2.9) 



We also have identified, by inspection, the sum over states and wave functions 
in Eq. (2.8) with the zero-frequency Lindhard-function Xo( r > r '; 0) of the non- 
interacting system. The introduction of an occupation number factor n(j) = 
1 — n(j) above is legitimate since all terms in the double sum where j is an 
occupied state cancel out due to the antisymmetry of the energy denominator. 
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Hence, the Newton- Raphson procedure suggests a density correction <5p( fc )(r) 
given by the linear integral equation (2.5), where F (r) is the functional 
(2.1) and 



e(r, r'; 0) = 5(r - r') - J d 3 r" Xo (v, r"; 0)y p _ h (r", r') (2.10) 
is the static dielectric function of a non-uniform electron gas [9] . 



3 Numerical Considerations 

A number of simplifications that are independent of the geometry can be 
applied to the algorithm to speed up the calculation; the justification for such 
simplifications is based on the observations that 

• The equilibration of the configuration is dominated by the Coulomb term 
Vc(r). 

• The dielectric function e(r, r'; 0) does not need to be evaluated very accu- 
rately for the iterations to converge, and 

• The ground state density is reasonably smooth. 

Generally, and in particular in the local density approximation, it is not a 
problem to calculate the variation <5Vn(r) / S p(r') . However, we found that the 
inclusion of the exchange-correlation term SV xc (r)/Sp(r') does not noticeably 
improve the convergence of the procedure: It is sufficient to use the Coulomb 
interaction 

U^')«j^. (3-1) 

This simplification allows us to use the same algorithm in microscopic calcula- 
tions [3] for which the computation of the variation of the exchange-correlation 
energy with respect to the density is more time consuming. Since the Coulomb 
term is vastly dominant, the non-local Fock term could also be included in 
the same manner and should not cause difficulties that are worse than the 
iterative solution. 

Since the kernel e(r, r'; 0) does not need to be very accurate for the convergence 
of the procedure, there is also no need to recalculate it after every iteration. 
This is especially true if one is already close to the ground state density. In 
fact, one can save the L — U decomposition of e(r,r';0), and needs to carry 
out only the back-substitution step during every iteration. Moreover, during 
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all those iterations when e(r, r';0) is not recalculated, Eq. (1.1) needs to be 
solved for the occupied states only. 

Finally, because the ground state density should be reasonably smooth, it is 
unnecessary to include very high lying single-particle states in the state sum 
(2.8). In fact, since the electron density falls off exponentially into the vacuum, 
it is in most cases sufficient to restrict the sum over particle states to those 
with negative energies. Even for very crude initial densities, we found that it 
is safe to include states with energies no higher than 0.2 Ry. 

With these considerations taken into account, the solution of Eq. (1.1) for any 
trial density p(r) is the most time consuming part of the algorithm. 

We conclude this section by remarking that the iterative solution conserves 
the particle number due to 



It is therefore an important precaution to start the iterations with an initial 
density p(r) that has the correct particle number. 



4 Application for spherical jellium clusters 

The feasibility of the calculation depends, of course, on the possibility to 
exploit the symmetries of the ground state and to reduce the three-dimensional 
eigenvalue problem (1.1) to a simpler task. If this can be accomplished, the 
calculation of the static dielectric function e(r, r'; 0) and the implementation 
of the Newton-Raphson iterations is equally feasible. 

We have implemented the algorithm for the symmetric slabs of jellium studied 
in Ref. [3], and for spherically symmetric jellium clusters which have in the 
past decade received much attention [10-12]. The convergence properties of 
the procedure are basically the same; we therefore describe the details for the 
cluster problem. 

In this geometry, the states 0i(r) are characterized by a radial quantum num- 
ber % and the angular momentum £, i.e. i = {i,£}. We consider closed-shell 
clusters, hence each state with angular momentum £ is 2{2£ + l)-fold degene- 
rate. Since the system is spherically symmetric, the Hartree- (or Kohn-Sham) 
equations decouple 




(3.2) 



<t>i (r) = <f> it (r)Y e , m (e,<l>) 



u it e (r) 



Y e , m (6,<j>). 



(4.1) 



r 
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Eq. (1.1) is a radial Schrodinger equation (note that we work in atomic units) 
" + + V H [p(r)] Ui , e (r) = e i/U u(r) (4.2) 

and the density is 

p(r) = i^(M)(2l + l)fo(r)| 2 . (4.3) 



The calculations follow exactly those of the derivation of the previous section; 
all angular integrations can be carried out due to the spherical symmetry. The 
angle-averaged Coulomb potential is 

Vc(r, r>) = f = ^ (4.4) 

J 4-7T |r — r I max(r, r'j 

and the angle-averaged static Lindhard function is 

dtt 



Xo(r,r';0) = J — Xo (r,r';0) 



£(2£ + l)n(i, i)n(j, e ^Ar)he(r)Mr')ht(r') ^ 



4vr 2 frr ' ' ' ' s j{ - z Ll 



In this geometry, the Newton-procedure becomes a matter of solving a one- 
dimensional linear equation. 

We demonstrate the working of our algorithm for the jellium model for Na 
clusters with 40 and 2018 electrons. To facilitate the comparison with earlier 
work, we use the local density approximation and the energy functional of Ref. 
[13]. The jellium background density p+(r) was taken to be a step function; 
as a worst-case scenario we assumed, as starting density, the same density 
profile for the electrons. During the iterations, the static susceptibility was 
re-calculated for the first two steps only and then kept fixed. 

Fig. 1 shows the density corrections Sp(r) during the first five iterations for 
the 2018 electron model; results for the average density correction as well as 
the kinetic, potential, and exchange-correlation energy are given in Table 1 
as a function of the iteration number. These results agree, within numerical 
accuracy, with those found in the literature [7,12]. Further documentation of 
the algorithm is provided in Fig. 2 where we show the difference between the 
energy obtained in the k-th iteration and the converged result. Also shown is 
the rms error in the density as function of the iteration number. Obviously, 
the rate of convergence is practically independent of the cluster size, and 
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comparable for different quantities of physical interest. We re-iterate that the 
replacement of the Coulomb potential by the full particle-hole interaction 
did not improve the convergence noticeably, updating the static susceptibility 
after every step reduced the number of iterations needed to obtain the accuracy 
shown in Fig. 1 to 16. 



5 Summary 

We have described in this note a rapidly converging algorithm for solving the 
Kohn-Sham and similar equations. While we have focused here on electronic 
systems, the algorithm is not restricted to those. Practically the same pro- 
cedure has been implemented in our work on bosonic liquid films [14] and 
droplets [15], the only difference being that only one state is occupied in 
that case. In self-bound systems, the "particle-hole interaction" cannot be 
replaced by the Coulomb potential, but otherwise the convergence properties 
of the Newton-Raphson algorithm are the same. 

For the implementation of our method, we have relied heavily on "canned" 
software, specifically BLAS and LAPACK routines for solving the band-sym- 
metric eigenvalue problems (1.1), (4.2) and for computing the L-U decompo- 
sition of the static dielectric function. The most time-consuming part of the 
algorithm is the solution of the eigenvalue equations (1.1). If needed, much 
computation time can be saved by choosing specialized algorithms, such as 
the "imaginary timestep method [16]", inverse iterations [17] or the Lanczos 
algorithm [18] that compute only the lowest lying states and work particularly 
well when good estimates for these states are known. 

We hasten to point out that none of the measures to speed up the computation 
is essential for the implementation of our algorithm for the simple systems 
considered here that can be reduced to one-dimensional problems. Even when 
the static dielectric function is updated every time, each iteration is a matter of 
a few seconds on a reasonable personal computer. Efficiency considerations will 
become truly useful only for higher-dimensional problems like non-spherical 
clusters. 

Finally, we emphasize again that a central quantity of our algorithm is the 
static dielectric function e(r, r'; 0) which is an independently interesting phys- 
ical quantity that one might wish to compute anyway. Our algorithm also 
highlights the reason for the delicate convergence of the normal iteration 
method: the dielectric function is vastly dominated by the potential energy 
term / d 3 r"xo(r, r'; 0)V r p _h(r", r') whereas this term is neglected by the con- 
ventional iterative procedure. 
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Fig. 1. The figure shows the density corrections 5p(r)/po for the electron density 
of a jellium cluster with r s = 4.0 and N = 2018 electrons during the first five 
iterations. The density is normalized to the density po of the jellium background, 
and the starting density was the same as the background density. The jellium edge 
is at r = 50.55 a.u.. 
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Fig. 2. The figure shows the relative energy corrections AE/Eoo = /E^ — 1 

where E^ is the converged result, and E^> is the energy obtained in k-th iteration, 
and the rms error e = y/J2i \^P( r i)\' 2 /(^Po) of the density. Results are shown for 
two jellium cluster with r s = 4.0, for N = 2018 (solid lines) and N = 40 (dashed 
lines). In both cases, the starting density was the same as the background density. 
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Table 1 

Convergence behavior of our iterations for a spherical jellium cluster with r s = 4 
and N = 2018 electrons. All energies are given in Rydberg units per electron, e is 
the rms error e = y/J2i \$P( r i)\ 2 / {Npo), the are the mesh points, iV is the number 
of mesh points, and po the jellium density. 
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